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ABSTRACT 

Context. Oscillations of stellar p modes, excited by turbulent convection, are investigated. In the uppermost part of the solar convection zone, 
radiative cooling is responsible for the formation of turbulent plumes, hence the medium is modelled with downdrafts and updrafts. 
Aims. We take into account the asymmetry of the up- and downflows created by turbulent plumes through an adapted closure model. In a 
companion paper, we apply it to the formalism of excitation of solar p modes developed by Samadi & Goupil (2001). 

Methods. Using results from 3D numerical simulations of the uppermost part of the solar convection zone, we show that the two-scale mass-flux 
model (TFM) is valid only for quasi-laminar or highly skewed flows (Gryanik & Hartmann 2002) and does not reproduce turbulent properties 
of the medium such as velocity-correlation products. We build a generalized two-scale mass-flux Model (GTFM) model that takes both the 
skew introduced by the presence of two flows and the effects of turbulence in each flow into account. In order to apply the GTFM to the solar 
case, we introduce the plume dynamics as modelled by Rieutord & Zahn (1995) and construct a closure model with plumes (CMP). 
Results. The CMP enables expressing the third- and fourth-order correlation products in terms of second-order ones. When compared with 
3D simulation results, the CMP improves the agreement for the fourth-order moments by a factor of two approximately compared with the use 
of the quasi-normal approximation or a skewness computed with the classical TFM. 

Conclusions. The asymmetry of turbulent convection in the solar case has an important impact on the vertical- velocity fourth-order moment, 
which has to be accounted for by models. The CMP is a significant improvement and is expected to improve the modelling of solar p-mode 
excitation. 

Key words, convection - turbulence - sun: oscillations 



1. Introduction 

In the uppermost part of the solar convective zone, turbulent 
entropy fluctuations and motions of eddies drive acoustic os- 
cillations. 3D numerical simulations of the stellar turbulent 
outer layers have been used to compute the exci t ation rates 



of solar-like oscillation modes Nordl und & Ste in (2001). As 
an alternative approach, semi-analytical modelling can pro- 
vide an understanding of the physical processes involved in 
the excitation of p modes: in this case, it is indeed rather 
easy to isolate the different physical mechanisms at work in 
the excitation process and to assess their effects. Various semi- 
analytical approaches have been developed by sev e ral authors 
dGoldreich & Keelevl Il977t iGoldreich et al.1 1 19941 iBalmforthl 
19921: ISamadi&Goupilll200lb : thev differ from each other 
by the nature of the assumed excitation sources, by the as- 
sumed simplifications and approximations, and also by the 
way the turbulent convection is described (see the review by 
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Stein et all 20041). Among the di fferent theoretical approaches, 
that of ISamadi & Goupi lf (l200ll) includes a detailed treatment 
of turbulent convection, which enables us to investigate differ- 
ent assumptio ns about turbulen t convection in the outer lay- 
ers of stars (Sama di et al. I l2005l) . In this approach, the analyti- 
cal expression for the acoustic power supplied to the p modes 
involves fourth-order correlation functions of the turbulent 
Reynolds stress and the entropy source term, which for the sake 
of simplicity are expressed in terms of second-order moments 
by means of a closure model. 

The most commonly used closure model at the level 
of fourth-order moments (FOM) is the Quasi-Normal 
Approximation (QNA), which is valid for a Gaussian proba- 
bility distribu tion function (se e Lesieurj[l997 ) and was first in- 
troduced by iMillionshchikov ( 194 lh . The QNA is r ather sim- 
ple and convenient to implement. However. lOgural dl963) has 
shown that such a closure could lead to part of the kinetic en- 
ergy spectrum becoming negative. In th is paper, we confirm 
the results of iKupka & Robinson! d2007h (hereafter KR2006), 
namely that this approximation indeed provides a poor descrip- 
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tion of the physical processes involved in solar turbulent con- 
vection. 

Mass flux models (e.g., Randall et al.1 [1992, 



Abdella & McFarland 11997b explicitly take the effects of 
updrafts and downdrafts on the correlation products into 
account. The presence of two well-defined flow directions 
then introduces an additional contribution when averaging 
the fluctuating quantities, since averages of fluctuating quan- 
tities over each individual flow differ from averages over 
the total flow. For applications in atmospheric sciences, the 
mass -flux model for convection has recently been improved 
by iGryanik & Hartmannl d2002l hereafter GH2002). Their 
motivation has been to account for the fact that horizontal 
scales of temperature and velocity fluctuations are different 
(hence their improvements lead to a 'two-scale mass-flux 
model' (TFM)) as well as to understand and measure the 
effects of the skewness of their distribution. According to 
GH2002, mass-flux models, which also include the TFM, 
underestimate the FOM by as much as 70%. Therefore, 
such models clearly miss some important phys i cal eff ects 
present in convective flows. IGryanik & Hartmann (2002) and 
Gryanik et alJ (120051) studied the asymptotic limits of TFM 
which led the authors to propose an interpolation between the 
QNA and the limit of large skewness provided by the TFM. 
This new parametrization permits a much better description 
of the FOM for convection in the atmosphere of the Earth 
(GH2002). We show that for their parametrization to be appli- 
cable to the case of solar convection, a more realistic estimate 
for the skewnesses of velocity and temperature fluctuations is 
required than that provided by the TFM itself (Sect. [2]). 

The parametrization of GH2002 requires the knowledge of 
the skewnesses and second-order moments to compute FOM. 
These have to be provided either by measurements, by another 
model, or by numerical simulations. In the present paper we 
do not aim to construct a complete model to compute these 
quant ities, which is the goal of the Reynolds st ress approach 
(e.g.. ICanutolll992t ICanuto & Dubovikov|[l998l) . Rather, we 
aim to analyze the shortcomings of the TFM and suggest im- 
provements using numerical simulations of solar convection as 
a guideline. The conclusions drawn from this analysis are used 
to derive a model for fourth-order moments in terms of second- 
order moments that can be used in computations of solar p- 
mode excitation rates. 

To proceed with the latter, we developed a formulation of 
the TFM that takes the effects of turbulence in each flow into 
account. This generalized TFM model (hereafter GTFM) is 
useful for both the superadiabatic and adiabatic outer solar lay- 
ers. This formulation can actually be applied in other contexts 
than just the excitation of solar p modes as long as the convec- 
tive system is composed of two flows. 

The GTFM is more general and realistic than the TFM, but 
it requires the knowledge of additional properties of both the 
turbulent upwards and downwards flows. We choose to deter- 
mine these properties by means of a plume model. Turbulent 
plumes are created at the upper boundary of the convection 
zone, where radiative cooling becomes dominant and where 
the flow reaches the stable atmosphere. In this region the 
updrafts become cooler and stop their ascent. This cooler 



flow is more dense than its en vironment and it trigger s the 
formation of turbulent plumes ( Stein & Nordlu nd 1998|). As 
shown by iRieutord & Zahn (119951) 7 these structures drive the 
dynamics of the flow; hence, to construct a closure model, 
we st udy the plume dynamics developed by IRieutord & Zahn 
( 19951) (hereafter RZ95). This makes it possible to build a 
closure model with plumes (CMP), which is valid in the so- 
lar quasi-adiabatic con vective region. In a companion paper 
(Belkacem et al. 2006, hereafter Paper II), we generalize this 
one-point correlation model to a two-points correlation model 
and calculate the power injected into solar p modes. 

The paper is organised as follows: Sect. 2 introduces the 
TFM. Its validity is then tested with a 3D numerical simu- 
lation of the uppermost part of the solar convection region. 
In Sect. 3, we extend the TFM formulation (GTFM) in order 
to take into account turbulent properties of both upward and 
downward flows. We next investigate the asymptotic limits of 
the GTFM. In Sect. 4, we construct the CMP with the help 
of the RZ95 plume model. We test the validity of this model 
with results from the 3D simulation and show that the use of 
the plume model limits the validity of the CMP to the quasi- 
adiabatic zone. The CMP is then used to obtain analytical ex- 
pressions for the third and fourth moments. Section 5 is dedi- 
cated to discussions and conclusions. 

2. The two-scale mass-flux model 

2.1. The model 

The TFM considers a convective medium composed of upward 
and downward flows that are horizontally averaged. The pres- 
ence of two flows introduces the possibility of a non-zero skew- 
ness for the moments of turbulent quantities when averages are 
done globally over the whole system. The TFM was developed 
in order to take into account this non-zero skewness. 

Any averaged turbulent quantity <p can be split into two 
parts, one associated with the updrafts and the other with the 
downdrafts: 



< (p > = a < <p >„ + (1 - a) < (p > d 



(1) 



where <> denotes ensemble spatial (in the horizontal plane) 
and time averages. < <p >» an d < <P >d are the averages for 
the upflow and downflow, respectively, a and 1 - a are the 
mean fraction al area occupied by the updrafts and downdrafts , 
respectively dRandall et al. 1 19921: IGryanik & Hartmannl 12002: 
Canuto & Dubovikovll998l) . 



Fluctuating quantities defined as <p' — <p- < (p > can be 
rigourously written as: < <p>' >- a < <f>' > u + (1 — a) < <p' >d, 
where the subscripts u and d are meant for upflow and down- 
flow, respectively. For vertical velocity fluctuations w', one 
then writes: 



< W > = a <w' > u +(1 - a) <w' >d 



(2) 



GH2002 propose to make the same decomposition for temper- 
ature fluctuations (6'); thus, hot and cold regions are considered 
separately. This step was motivated by the observation that for 
the case of atmospheric boundary layer convection the charac- 
teristic horizontal scales of velocity and temperature flucuta- 
tions are different from each other and by the fact that the plain 
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Fig. 1. On the top, the superadiabatic gradient (V - V a d) is plot- 
ted versus the depth (z). The reference depth (z = Mm) corre- 
sponds to the photosphere. At the bottom, the mean fractional 
area of the upfiow (a) and the warm drafts (b) are given. To cal- 
culate these quantities the upfiow and downflow are separated 
using the sign of w' as a criterion. The same is done for the 
warm and cold drafts. 



mass flux average Eq. (Q]l violates certain symmetries between 
velocity and temperature fiucutations. Indeed, hot and cold re- 
gions do not necessarily coincide with updrafts and downdrafts, 
respectively. Hence, a second quantity (b), the mean fractional 
area occupied by warm drafts, is introduced, and in most cases, 
a + b (thus the name TFM). Then, 



< ff >= b < ff > h +(1 - b)<& > c 



(3) 



Furthermore, the TFM defines the velocity fluctuations in- 
side the upfiow (w' u ) and downflow (w' d ), respectively, as: 

w' u — w u — < w > and w' d — wj— < w > . (4) 

Similarly, for the temperature fluctuations inside hot (ff h ) and 
cold (9' d ) regions, respectively, one has 

0' h = 6 h - <6> and ff c = G c - < 6 > . (5) 

The quantities w u , Wd, 8h, and 6 C are the averages of veloc- 
ity and temperature, respectively, over all updrafts (w„), down- 
drafts (wd), hot (Oh) drafts, and cold (8 C ) drafts. Clearly, aver- 
ages of the four fluctuating quantities in Eqs. © and (0 do not 



Fig. 2. The skewnesses S K (on the top) and Sg (on the bottom) 
are plotted versus the depth (z). Solid lines represent direct cal- 
culation from the 3D numerical simulation (Eq. Q) and dashed 
lines represent the skewnesses calculated using the TFM model 
(Eq. (0). 



vanish because the average of a quantity over the whole flow 
differs from the average over one single (up or down, hot or 
cold) draft. 

It is expected that the differences between the updrafts and 
downdrafts lead to a probability distribution function (PDF) 
that is no longer symmetric with respect to vanishing veloci- 
ties and temperature differences. As the third-order moments 
(< w' 3 > and < 6' 3 >) vanish when the PDF is symmetric, their 
values provide a measure for the deviation from a symmetric 
PDF. The skewnesses are defined as: 

^ ,,/3 ^ ^ a'3 ^ 



S w — 



and Sg — 



,3/2 



(6) 



respectively, for the vertical velocity and temperature fluc- 
tuations. In order to compute expressions for higher order 
moments in terms of velocity and temperature fluctuat ions, 
Eqs. © and ©, GH2002 followed iRandall et all dl992l) . us- 
ing an additional simplifying approximation, i.e., 

<4>" > ~ < >" , (7) 

9', }. This approximation neglects the con- 



VI' 



where <p 

tributions of fiucutations within the up- and downdrafts and 



w,d' h.c 
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differences in temperature and velocity between the individual 
drafts. 

Given this approximation and the known second-order mo- 
ments, the TFM provides third-order moments as follows (see 
GH2002): 



< w ,2 e' > = s w < w' 2 > 1/2 < w'e' > 

< w'e' 2 > = s 9 <o j2 > 1/2 < w'ff > 



(8) 



and FOMs as: 



<w' 4 > 



(1+S 2 ) <w' 2 > 2 



9' 4 > = (1+S 2 ) < 



f 2 >2 



(9) 



< w'V > = (1 + S 2 J < w' 2 > < w'. 

< w 'e' 3 > = (i + sh < e' 2 > < w'e 



> 
> 



The skewnesses S w and Sg (Eq. (|6]i) are related to a and b 
through 



S w 



1 -2a 



and S 



I -2b 



o 



(10) 



V«(l -a) ^b(l-b) 

(GH2002, see also lRandall et al.lll992l for the case of S w ). 

In the following we consider only vertical-velocity mo 
ments. Assuming S w — S$ — in Eq. (0 gives: 



>4 >2 2 

<W > = < W > 



(ID 



Such a result is not consistent with a quasi-normal (Gaussian) 
PDF. In deed, when w' follows a normal distribution, then 
dLesieujll997b : 



and 



<w' 4 > 



3 < w' 2 > 2 



(12) 



GH2002 found that the two-scale mass-flux average, Eqs. dS),- 
,([T0ll, underestimates both skewness and fourth-order moments 
as measured by aircraft data for planetary boundary layer con- 
vection (see their Figs. 4 and 7). To account for the omitted 
contributions from fluctuations within and between the up- and 
downdrafts, they suggested generalizing the TFM by build- 
ing the fourth-order moments as an interpolation between two 
asymptotic regimes: 

- S w — 0, assuming the quasi-normal approximation (QNA) 
limit that is valid for a Gaussian PDF, and 

- S w » 1, the large skewness limit (GH2002). 

GH2002 hence proposed: 



<w' 4 > 



< w'V > 



< w'0' 3 > 



3 (1 + <W 2 > 2 

3(1 + ^S 2 ) <0' 2 > 2 

3(1 + is 2 ) <w' 2 > <w'6'> 

3 (i + -sj) <e' 2 > < w'e" > 



(13) 



Corresponding expressions for other FOMs (< w 6 > 
and those includin g horizontal velocities) can be found in 
Grvanik et al.1 d2005l hereafter GH2005). 




z (Mm) 

Fig. 3. Fourth-order moment (< w' 4 >) as a function of depth 
(z) normalized to the FOM, as calculated directly from the sim- 
ulation. The solid line denotes the moment calculated using 
Eq. ( fT3l with S w taken directly from the simulation; the dashed 
line shows the result if S w is instead taken from Eq. (fTOb . as 
in the TFM case; and the dotted line is the QNA (Eq. (|26||). 
Equations. dT3~b and d26b involve second-order moments that 
are computed using the numerical simulation. 



2.2. Validation with a 3D numerical simulation of the 
solar external layers 

We consider the uppermost part of the solar turbulent convec- 
tion. Turbulent plumes are known to exist wi thin this region 



(Catt aneo et all 1 1 99 U IStein & Nordlund 1998). Here, we test 



the validity of the TFM using 3D numerical simulations of 
these upper solar layers. The geometry is plane-parallel with a 
physical size of 6 Mmx6 Mmx3 Mm. The upper boundary cor- 
responds to a convectively stable atmosphere and the lower one 
to a quasi-adiabatic convection zone. The 3D simulations used 
in this work were obtained with Stein & Nordlund's 3D nu- 



merical code ( Stein & Nordlundl ll998). Two simulations with 
different spatial grids were considered: 253 x 253 x 163 and 
125 x 125 x 82. 

Averages and moments of the velocity and temperature 
fluctuations were computed in a two-stage process: 
a is given as the number of grid points per layer with upwards 
directed vertical velocity divided by the total number of points 
in that layer. The instantaneous value of b is obtained in a 
similar manner, comparing the temperature at a given point in 
a layer with its horizontal average. Moments related to updrafts 
were obtained from horizontal averaging, using only those 
grid points at which vertical velocity was directed upwards 
at the given instant in time, and likewise, quantities related 
to downdrafts were obtained from horizontal averaging using 
only those grid points at which vertical velocity was directed 
downwards. In a second step, time averages were performed 
over a sufficiently long period of time such that averages no 
longer depended on the integration time beyond a few percent. 

— Calculation of the skewnesses: Computations of the 
mean fractional area of the upflow (a) and downflow (1 - a), as 
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well those of the warm (b) and cold (1 — b) drafts from the nu- 
merical 3D simulations (Fig.[TJ, show that the upper part of the 
solar convection zone can be divided into three parts: the stable 
atmosphere, the superadiabatic zone, and the quasi-adiabatic 
zone. In the convectively stable atmosphere (z < Mm, where 
z — is approximately at the bottom of the photosphere and 
z = -0.5 Mm is the uppermost boundary of the simulation), 
there are no asymmetric motions. In the superadiabatic zone 
(0 < z < 0.5 Mm), from the top downwards, the depar- 
ture from symmetry for the flows strongly increases (Fig. [T), 
and the skewnesses, S w and Sg, significantly differ from zero 
(Fig. |2). Hence, one must expect a non-negligible departure 
from the QNA, which is explained by radiative cooling creat- 
ing turbulent plumes. In the quasi-adiabatic zone, plumes have 
already been formed and no additional asymmetry is there- 
fore created. Hence, the asymmetry remains large and constant 
(a =s b x 0.7) and the skewnesses show a constant departure 
from S w = S e = 0. 

The last two regions are of interest in this work because 
both show a departure from the quasi-normal PDF in terms 
of fluctuating vertical velocity and temperature. The com- 
parison of the above numerical results with the results from 
the classical TFM (Eq. ©) and the TFM model (Fig. 0) 
shows that Eq. dTOb fails to reproduce the behaviour of the 
skewnesses from the 3D s imulation (as was also found by 
Grvanik & Hartmannl 12002 for convection in the atmosphere 
of the Earth, see their Fig. 4). 

— Detailed comparison of a fourth-order moment: The 
GH2002 interpolation relation Eq. ( fT3l combined with the 
TFM relation for skewness, Eq. dTOb . shows only a slight im- 
provement of the QNA description for the FOM < w' A >, when 
compared to the numerical result (Fig. [3]). 

To conclude, it seems that a physical process is missing 
in the quasi-adiabatic convective zone. To explain such a dis- 
agreement between the numerical results and the TFM, we 
must come back to its main approximation (see Eq. (0). For 
n — 2, Eq. (0 yields: 



< w' 2 > - < w > 2 











Hence, the TFM assumes that the variances of the fluctuations 
of vertical velocity and temperature within and among individ- 
ual drafts vanish, and the detailed turbulent nature of the flows 
themselves does not have to be taken into account. In order to 
compensate for the shortcoming of Eq. © and thus the conse- 
quences of the approx i matio n Eq. 01 on the model predictions, 
Grvanik & Hartmannl d2002l) proposed a more general interpo- 
lation relation (Eq. ( fT3l ) that uses Eq. (0 only for one of two 
asymptotic limits. 

As seen above, Eq. ( TTOb fails to describe the numerical re- 
sults. The question therefore is whether the interpolated rela- 
tion (Eq. ( fT3l ) is still valid, provided a correct value for the 
skewness is used. Hence, we assess the validity of Eq. (fT3l ) 
by inserting the value of S w directly given by the 3D numer- 
ical result. The result is shown in Fig. [3] as well. This is the 
model that IGryanik & Hartmannl (120021) proposed to be used 
instead of the TFM itself and its associated relation for the 
skewnesses, Eq. (TTOb . We obtain an accurate description of the 
FOM < w' 4 > in the quasi-adiabatic region, but not in the 



superadiabatic zone, where the interpolated relation does not 
seem well adapted (cf. KR2006 for a more detailed discussion). 

3. The generalized two-scale mass-flux model 

3.1. Theoretical formulation 

Here we remove the approximation of Eq. (0 and instead con- 
sider the exact expression: 



< w' n >= a < w' n >„ + (1 - a) < w' n > d 



(15) 



Our main idea is to separate the effect of the skewness in- 
troduced by the presence of two flows from the effect of 
the turbulence that occurs in each individual flow. We note 



that in a geophysica l context Sie besma & C uiipers (1995) and 
Petersen et al.1 (1 19991) studied the transport properties of classi- 
cal mass-flux models that also involved a separation of large- 
scale and turbulent components. Here, we start from th e more 
rece nt viewpoint of the TFM by Grvanik & Hartmannl d2002b 
and IGryanik et al.1 d2005l) . which takes into account that up- 
drafts and downdrafts are not strictly correlated with hot and 
cold drafts, respectively. As a first step we define the intrinsic 
fluctuations within one of the flows as: 



Wj = Wj- < w >j 



(16) 



where j = {u, d). They are fluctuations with vanishing averages. 
To express W- in terms of w'j (Eq. (0), we write: 



Wj — Wj + < W >j — < w > . 



(17) 



Applying the decomposition of Eq. (0 to < w > in the above 
expression yields: 



w u - w u + (1 - a) 6w 
w' d — w' d - a 6w 

with 



(14) 6w -< w >„ - < w > d - | < w > u | + | < w >d 



(18) 



(19) 



because < w > u > and < w >d < 0. 

Inserting Eq. (fT8l into Eq. <TT3T > for n — 2, 3 yields: 



< w > 



a(l - a) Sw 2 



+ a < w > u + (1 - a) < w >d 

< w' 3 > = fl(l -a)(l -2a)6w i 

+ a <w' 3 > u + (1 - a) < w' 3 > d 



(20) 



+ 3a(l - a) 



< W' 2 > u - < W' 2 >d 



5w . 



(21) 



The third-order moment (Eq. (f2Tli). which is related to the 
skewness (see Eq. (O), is composed of four contributions: 

- the first term i s the expression derived by 



Gryanik & Hartmann (2002). It is a measure of the 



skewness introduced by the presence of two flows. 
- the second and third terms represent the asymmetry of the 
PDF within each flow induced by turbulence. 
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- the fourth term measures the difference of the fluctuating 
velocity dispersion. Hence, if one of them is larger than the 
other, the PDF becomes asymmetric. 

The description of the turbulence in individual flows that has 
been neglected in the TFM is included in the present formula- 
tion through the last three terms in Eq. (f2TJ. 

We next focus on the fourth-order moment < w' 4 >, which 
is of interest in the context of stochastic excitation of solar 
p modes (see Paper II). Then setting n — 4 in Eq. ( fT5l l. we 
have: 



5w 



< w' 4 > = a{\ - a){\ -3a + 3a 2 ) 5w A 

+ 6a(l - a)((l - a) < w' 2 > u +a < w' 2 > d \ 

+ 4a(l - a)^ < w' 3 > u - < w' 3 >d jSw 

+ a < w' 4 > u + (1 - a) < w' 4 > d . (22) 

We stress that the TFM is recovered from the present gener- 
alized formulation when proper fluctuations (i.e., turbulence) 
within and among the individual drafts are neglected, i.e., 

< w' n >= 0. 

The same decomposition can be performed in terms of tem- 
perature fluctuations. As the calculation is symmetrical in w', a 
and ff , b, we hence have: 

<g' 2 > = b(l -b)S6 2 

+ b<e' 2 > h +(i -b)< ~e' 2 > c 



9' 3 > = b(l -b){\ -2b)50 3 

2/3 



+ b < S' 3 > h + (1 -b) < 6' 3 > 



3b(l-b) 



< e' 2 > h - < y 2 > c 



36 



9' 4 > = b(l - b)(l - 3b + 3b 2 ) , 



+ 6^(1 - b)((l -b)< S' 3 > h +b < S' 2 > c ) S9 2 
+ 4/7(1 - b)[ < ~0' 3 > h - < Q' 3 >c ) 68 



+ b<0' 4 > h +(1-/7) < <9' 4 > c 



S'4 



(23) 



The next step consists of the derivation of the cross terms < 
w'ff >, < w' 2 6' 2 >, < w' 2 ff > and < w'6' 2 >; it is convenient 
to define the coefficients a uc so as to take into account the 
four types of flow (see also GH2005): 

- warm updraft, a u i, 

- cold updraft, a uc - a- a u h 

- warm downdraft, a^h = b — a u h 

- cold downdraft, aj c = 1 - b - a uc 

Expressions for the third and fourth cross-correlation moments 
are given in AppendixlAl 

The generalized TFM has the advantage of isolating the 
skewness introduced by the two flows (as measured by S „, and 
Sg in Eq. ( [Tol l) from the effects of turbulence in each of the 
flows (as measured for instance by the two terms w'^ and w' 2 ). 
The GTFM allows us to take the effects of turbulence into ac- 
count. We note that a small value of the kurtosis can occur only 
if proper fluctuations lead to negligibly small deviations from 



the root mean square average. Such a flow pattern consisting 
of clearly defined up- and downflows as well as hot and cold 
areas with a kurtosis K w > 1 can be considered as representing 
a quasi-laminar state . We stress that for the quasi-laminar case, 
Eq. (0 remains exact; thus the kurtosis becomes: 



K w = 



< w' 4 > 

< w' 2 > 2 



= (1 +S 2 V ) with5 w = 



1 -2a 



(24) 



Vfl(l - a) 

For a = 0.5, one obtains K w = 1, which is far from the value 
for a Gaussian PDF (K w = 3). To take into account turbu- 
lence within the up- and downdrafts, one can use Eq. ( fT3l (see 
Sect. 12.2b with the skewness S w =< w' 3 > / < w' 2 > 3 ^ 2 from 
the GTFM. In this case we obtain: 

K w = 3(l + ^S 2 W ). (25) 

This implies that a (moderately small) non-vanishing skewness 
will make the value of K w closer to three than in the quasi- 
laminar case. In the solar case, in the quasi-isentropic zone 
S 1 ], ~ 4 (Pig. [2j, hence K w « 3 + 4/3. In the physical picture 
underlying Eq. (fT3l . turbulence prevents the PDF from being 
too far from a Gaussian one (K w — > 3). 

We notice that one important source of turbulence that can 
be considered responsible for at least part of the fluctuations in 
a draft — in addition to those created by the radiative processes 
on top of the convection zone — is related to shearing stresses 
between the up- and downdrafts. However, the investigation of 
the sources of turbulence is beyond the scope of the present 
work. Those mechanisms certainly play an important role in 
both the small scale velocity and the thermal fluctuations. Their 
study is definitely desirable in the future. One should also note 
that the splitting approach of the GTFM is valid and can be 
used for any convective system, provided that it is composed 
of two flows. As it is unclosed, it must be seen as a good basis 
for building a closure model. 

3.2. Asymptotic limits 

In the following, we study the asymptotic limits of the GTFM, 
focusing on the fourth-order moment < w' 4 >. The stan- 
dard mass flux model is easily recovered when setting the 
proper moments to zero: < w' >= in Eqs. ( 1201 .— . 1(221 . 
The same holds for the TFM, Eqs. (H),-,®, which is recov- 
ered, if in addition < §'" >- i n Eq. ( l23b (cf. Eqs. (7) and 
(8) in iGrvanik & Hartmannll2002h . We now turn to the QNA 
limit and the limit for large skewness, which are more inter- 



they c 
[al] d2 



esting as they are u sed by IGrvanik & Hartman n (2002) and 



Grvanik et alj (120051) in order to corroborate the interpolation 



formula Eq. ( fl3b . 

3.2.1 . The quasi-normal limit 

To obtain the QNA (Eq. (fl2l) ). it is necessary that S w = 0, but 
it is not sufficient. In fact, a vanishing skewness only shows 
that the PDF is symmetric, but not that the PDF is Gaussian. 
Further conditions are necessary: 

- the moments must have zero mean, which implies | < w >,, 
| = | < w >d I = from Eq. © and Eq. ( fT8l : 
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- for the QNA to apply to the whole system, one must assume 
that the QNA is valid for each flow; 

- we must also assume that a = 0.5; 

- the turbulent pressure must be the same in the upfiow and 
downflow. Otherwise the skewness (S w ) is different from 
zero, according to Eq. (l2TT i. and the consequence is an 
asymmetric PDF, which is not consistent with the quasi- 
normal assumption. So the condition < vv' 2 >«=< vv' 2 >d is 
required. 

Then starting with Eq. (l22l . we find: 

' 4 > = a < vv' 4 > u + (1 - a) < vv' 4 > d 



< w 



— < vv 

2 



,~>2 



H — < W 

2 



r'2 



and finally 



< w > - 



3 < w' 2 > 2 



(26) 



which is the expression for the fourth order moment in the 
QNA. Note that the TFM (Sect. is unable to properly re- 
cover the QNA. Within the GTFM the QNA results from two 
terms, < vv' 4 >„ and < vv' 4 >d, which are related to the intrin- 
sic turbulence in each flow, but these are neglected in the TFM. 
This example also demonstrates that for a convective flow, the 
deviation of a PDF from a Gaussian one cannot be modelled by 
the TFM without further modifications of that model (even if 
a = 0.5). 

3.2.2. The large skewness limit 

Grvanik et alj d2005l) have shown that the TFM must be recov- 
ered when considering a convective system with large skew- 
ness. Then, for S w » 1, the expression for < w' 4 > in Eq. ( fT3l 
becomes: 




z (Mm) 



Fig. 4. Second-order moment of the upfiow over that of the 
downflow (< vv' 2 > u I < vv' 2 ><*) as a function of depth, cal- 
culated directly from the simulation. Upfiow and downflow are 
determined according to the sign of vv' . 



interpolated expressions for the FOMs (Eq. ( fT3l ) are limit- 
ing cases for a flow that consists of a coherent part with two 
components (up- and downdrafts), which themselves are sub- 
ject to turbulence (cf. the discussion of the GH2002 model in 
KR2006). In Sect. I2.2l we have shown, using the 3D numerical 
simulation, that this interpolation is valid provided the skew- 
ness is taken directly from the 3D simulation. This property 
can be understood using the GTFM, as it permits us to ob- 
tain the different ing r edient s of the interpolation formula of 
iGryanik & Hartmann (2002) from Eq. (l22l and the individual 
contributions to Eq. d22l . can be analyzed using numerical sim- 
ulations. 



<w' 4 >*Si <w' z > 2 



(27) 



The large skewness limit physically corresponds to either a « 1 
or a « 0. Indeed, it means that one of the two flows dominates 
over the other one in terms of mean fractional area in the hori- 
zontal plane. Thus, due to conservation of mass, the mean ver- 
tical velocity becomes large such that 6w » 1 m s _I in the 
solar case (see Sect.@] Eq. (HTt ). 

In Eq. d22l) . the term proportional to 6w 4 , which measures 
the effects introduced by an asymmetric flow, dominates and 
leads to the TFM expression for the fourth-order moment < 



w' 4 >: 



/4 

< vv >= a 



(1 -a)(a 3 +(1 -a) 3 j(5w 4 



(28) 



Grvanik & Hartman n (2002) demonstrated that this expression 
leads to the relation: 



4. The closure model with plumes 



Section lZ2l confirmed the conclusion by KR200 6 that the inter- 
polate d relations in Eq. ( TT3l proposed bv lGrvanik & Hartmann 
(2002) could be adapted for the solar case provided that the 
skewnesses are appropriately calculated. Using the GTFM to 
model skewnesses, Eq. (f2TT > shows that the skewness S w , for 
instance, depends on six quantities: 5w, < vv' 3 > lu d, < vv' 2 > u ,d, 
and a. As shown below, some of the terms in S w turn out to 
be negligible in the quasi-adiabatic convective region because 
plumes are more tu rbulent in the downflow than in the upfiow 
(Stei n & Nordlundlll998l) . The remain ing dominant terms are 
modelled hereafter by a plume model jRieutord & Zahnll 19951) 
in the quasi-adiabatic convective region, where the CMP is 
valid. 



< w' 4 >= (1 + S 2 ) < w' 2 > 2 ~ s 2 < w' 2 > 2 for S » 1 , (29) 4.1. Turbulence in upflows and downflows 



where, as in Eq. ( [Tol l, S = (1 - 2d) / ~\Ja(l - a). The same would 
result if the exact function S w were taken in this limit instead 
of its approximation, Eq. ( [Tol l. 

Hence, the GTFM enables us to show th at the asymptotic 
limits used by IGryanik & Hartmann! (120021) to motivate the 



In Fig.|U we compare the second-order moments of both flows. 
These quantities are of the same order of magnitude in the up- 
per part, above the photosphere. From the photosphere, the ra- 
tio < vv' 2 >„ / < vv' 2 >d then sharply decreases, with increasing 
depth (z). Hence, contributions to the skewness (S w ), involving 
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z (Mm) 



z (Mm) 



Fig. 5. The terms 3a(l - a) < w' 2 >d Sw (solid line), a < vv' 3 >„ 
(dot-dot-dot-dashed line), and (1-a) < vv' 3 >d (dashed line) are 
plotted versus the depth (z). From Eq. (|2TI ). the dominant terms 
remains 3a(l - a) < w' 2 >d 6w. This justifies the assumptions 
that the terms involving third-order moments can be neglected 
in the quasi-adiabatic zone. 



Fig. 6. Mean velocity profile of the upflow (dashed line) and 
downflow (solid line) as a function of the depth. Note that the 
peak at z = 0.1 Mm corresponds to the maximum turbulent 
pressure. The use of power laws limits the validity of the CMP 
to the quasi-adiabatic zone, as is implied by the deviation of the 
profiles from power laws in the superadiabatic region. 



< w' 2 >„ (Eqs. (|20| > and Eq. ( T2TI )) can be neglected in com- 
parison with those involving < w' 2 >d in the quasi-adiabatic 
part of the convection zone. The third-order moments < vv' 3 >d 
and < vv' 3 >„ can also be discarded (see Fig. [5]) because their 
contributions are negligible. 

The skewness S w then becomes: 



S w = 



a(l — a) 

< w' 2 > 3 / 2 



2a)dw 2 -3 <w 



a > d ) 



6w . 



(30) 



where 6w is given by Eq. 119) . Hence, only < w' 2 > ( ; and 6w re- 
main to be modelled. Similarly, the 3D calculations show that 
the cool medium is more turbulent than the hot one and that 
third-order moments for the temperature fluctuations can be ne- 
glected. Then the expression for Sg becomes: 

S e = b{\ - b) <y2 1 >3/2 ((l - 2b)S0 2 - 3 < V 2 > c ) 69 , (31) 

where the quantities 69 =< 9 >/, - < 9 > c and < 9' 2 > c must 
be modelled. 

Note that in the QNA limit 6w = 0, so that for the expres- 
sion Eq. d30l >. S „ = 0, and according to Eq. ( fT3l . 

< w' 4 >= 3 < w' 2 > 2 . 

However, because we have assumed < w' 2 >„ << < w' 2 >d 
when deriving the expression S w , rigourously speaking, S w 
does not tend correctly to zero in the QNA limit. Such an ex- 
pression therefore cannot be used in the case of a near QNA 
regime. In our case, we have shown in Sect. l2.2l that the medium 
is far from the QNA limit in the quasi-adiabatic zone, and hence 
the expression Eq. ( f30b can be safely used. 

To proceed further, < w' 2 >d and < ff 2 > c are written in 
a more suitable form. We neglect < w' 2 >„ in Eq. d20b for 

< w' 2 >, and < 6 2 > h in Eq. (|23]) for < 9' 2 >. This yields: 



< w' 2 > - a(l - a) 5w 2 + (1 - a) < w' 2 > d 
<9' 2 > = b(l - b) S9 2 + (l-b)< 9' 2 > c . 



We then derive expressions for < w' 2 >d and < ff 2 > c in terms 
of < w' 2 >, 5w, and < 9' 2 >, 69, respectively (see Eq. (l32l). 
Inserting them into Eqs. (f30b and OTb . one then obtains: 

S w = l — — ■ a ((1 - fl )(l - 5a)6w 2 - 3 < w' 2 > )6w (33) 

< w' 2 > 3 ' 2 \ ) 

and 

S g = 1— b ((1 - b){\ - 5b)69 2 - 3 < 9' 2 > )&9 . (34) 

< 9' 2 > i/2 \ ! 

We assume that the second-order moments (< w' 2 > and 
< 9' 2 >) are known. In the present work, they are computed 
from the 3D numerical simulation. In principle, they could also 
be taken from a convection model such as the mixing-length 
theory. The last step then is to determine 6w and 69 (as well as 
a and b). As 6w is the difference between the mean velocities of 
upward and downward flows, it is possible to model it by means 
of a plume model. This approach is also used to determine 69. 



4.2. The plume model 
4.2.1. Determination of Sw 



We u se the model of plumes developed by iRieutord & Z ahn 
(119951) . The plume is considered in an axisymmetric geometry 
with a Gaussian horizontal profile for the vertical velocity (wv), 
the fluctuations of enthalpy (6h), and density (6p) such that 

w d (r,z) = V(z) exp(-r 2 /b 2 p ) , 

6p(r,z) = Ap(z) exp(-r 2 /b 2 p ) ,and 

6h(r,z) = M(z) exp(-r 2 /^) . (35) 

where b p (z) is the radius of the plume. We assume, as in RZ95, 
an isentropic and polytropic envelope structure, hence 



(32) 



p(z) = po(zlz Q ) q 
T(z) = To(z/z ), 



(36) 
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Table 1. Solar values of plume model parameters (from RZ95) 





* 0.1 


Pa 


190 kg m- 3 


/• 


3.9 10 26 W 


Zo 


« 2.10 8 m 


go 


270 m s~ 2 



where q is the polytropic coefficient, po and Tq are the density 
and temperature at depth z = Zo, and zo is the reference depth 
that corresponds to the base of the convective region. 

In Fig. [6] we show that the mean velocity of upflow and 
downflow in the quasi-adiabatic convection zone both obey a 
power law in (z/zo)' - We therefore assume a power law for the 
mean velocity of the downflow (i.e., the plumes). Then 



<w> d = w d o 



( £ )' 



(37) 



with 



w d0 



12F 



1/3 



(38) 



(RZ95), where r = {-q + l)/3, /3 = 3a/(q + 2), and a 



0.083 is the entrainment constant for a Gaussian profile (Turner 
1986). F is the convective energy flux and g is the gravitational 
acceleration. In Table [TJ we list solar values of the previously 
introduced parameters taken from RZ95. These values are used 
in the present paper except for F, which is taken from the 3D 
numerical simulation (as explained below). For a monoatomic 
perfect gas, one has q = 3/2, hence r = -1/6. However, our 
3D numerical simulations indicate a value of r closer to 0. The 
reason is likely that there is radiative cooling. Hence, y > F = 
Cp/cy, where y is the polytropic index (q = l/(y — 1)). 

Following iRieutord & Zahnl (11995b . we assume that all the 
convective energy flux is transported by the plume, thus 



F = L Q IN , 



(39) 



where N is the number of plumes in the shell at h — R Q — z- We 
find N ~ 6.10 6 from the 3D numerical simulation. To obtain 
such a result, one has to use the relation between a and N: 



a = Nnb 2 p /4nh 2 



(40) 



where (a) is mean fractional area of the upflow, h — R Q — z, and 
b p is the radius of a plume. (b p ) and a are taken from the 3D 
numerical simulation. We assume a = 0.7, as taken from Fig.Q] 
which shows that the mean fractional area a is roughly constant 
in the quasi-adiabatic convection zone. The plume radius, b p , 
is estimated at the top of the simulated box, which corresponds 
to the photosphere. 

At this stage, we have modeled the downdrafts, but not yet 
the updrafts. The 3D numerical simulations show that mean 
velocities of upflow and downflow obey the same power law 
(Fig. [6j. This can be explained as follows: from the conserva- 
tion of the mass flux one has 



< pw > — a < pw > u + (1 - a) < pw >d = . 



(41) 




z (Mm) 

Fig. 7. Fourth-order moment < w' 4 > as a function of depth 
z normalized to the FOM calculated directly from numerical 
simulations. The solid line shows < w' 4 > calculated using the 
CMP model, the dashed line is the moment as obtained from 
Eq. (O with Eq. dTOj for S w , and the dotted line is the QNA, 
Eq. (HQ. 



Fluctuating parts of densities in up and downflows are ne- 
glected such that p u ;»< p >„ and p^ «< p >d (see Fig. 2b 
ofRZ95).Thus, 



< w > u 



(l-a)<p > d 
a < p >„ 



< w >d 



(42) 



Then, assuming that < p >„= p u o(z/zo) q and < p >d- 
Pdo(z/zo) q obey the same power law as in Eq. 



< W > u 



(1 - a)pdo 
a Pm 



<w> d . 



(43) 



PuO and pdo are the values at the reference depth zo- We set 
a xs 0.7 (see Fig. [TJ, which is the value obtained in the quasi- 
adiabatic zone from the 3D numerical simulation. Assuming 
further that pdo/p u o ~ 1, one obtains 



I <w >„ I ~ 0.45 | < w > d | , 

which is approximately what is seen in Fig. [6] 



(44) 



4.2.2. Skewness S w and the fourth-order moment 

We use Eq. (1331 for the skewness with 

Sw = (< w u > - < w d >) w 1.45 w dQ . (45) 

The vertical depth of the computation box is narrow in com- 
parison with the reference depth zo, thus 6w varies only weakly 
with z- Hence, we assume r — in the solar case. The fourth- 
order moment < w' 4 > can then be computed by means of the 
interpolated relation Eq. ( TT3l . In Fig. [TJ we show the resulting 
< w' 4 >. The CMP clearly is an improvement compared to the 
QNA and the TFM expression for S w , Eq. ( TTOb combined with 
Eq. (fT3l , by at least a factor two in the quasi-adiabatic zone. 
The FOM in the superadiabatic zone is overestimated. Indeed, 
as mentioned above, the CMP is not able to describe such a 



10 



K. Belkacem et al.: A closure model with plumes 




z (Mm) 

Fig. 8. Fourth-order moment < ff A > as a function of depth 
Z normalized to the directly numerically calculated FOM. In 
solid lines the moment stems from < 9' 4 > calculated using the 
CMP model, the dashed line is the moment as obtained from 
Eq. ( f]~3T > and Eq. ( TTOt for S e and the dotted line is the QNA. 



zone mainly because the assumptions of Sect. 14. H are not valid. 
Note that it is possible to use the same procedure to compute 
any other third- and fourth-order moment. 

4.2.3. Determination of 69 

Similarly to the procedure in the previous section, we evalu- 
ate < 9' 2 > with the help of Eqs. (f32| and (O. We therefore 
need to determine 59. The temperature profile is more sensi- 
tive to departure from adiabaticity than the velocity profile. It 
is therefore not suitable to assume an isentropic envelope. Such 
an approximation can still be used in the downflow, but not for 
the upflow, which is far from being adiabatic due to radiative 
cooling. Then, for the sake of simplicity, we assume a power 
law to obtain 59: 



59 * 69 (-)" 



(46) 



For z > 1 Mm in the simulated box (z — Mm denotes the 
photosphere), one derives m = -1.5, 69q « 170 K from the 
3D numerical simulation. Using the power law (Eq. d46l l) with 
m = -1 .5, the skewness Se can be calculated using Eq. 1341 . In 
Fig. [8] we present the fourth-order moment < 9' 4 > computed 
using the CMP, and as expected, the description of the FOM is 
improved. In the deeper part of the convection zone (i.e., the 
adiabatic region), 59 is easier to model because Eq. d36l l can be 
used and the difference 59 becomes a power law. From Eqs. (JHJ, 
( fT3l ). d33l . and ( 134b all the third- and fourth-order moments can 
be modelled with the CMP. 

4.2.4. Summary: the CMP in a nutshell 

In practice, one uses the CMP to compute < w' 4 > by means 
of the interpolation formula Eq. ( fT3b . where the second-order 
moment < w' 2 > is supposed to be known and where the skew- 
ness S w is computed from Eq. d33l . In the latter expression, 5w 



is determined using the plume model through Eq. d43T > and us- 
ing Eqs. d3~7b to (l40l > with appropriate values of parameters for 
the case studied (in the present paper we used the values from 
Table [1] which are suitable for the solar case). Here, a(z), N, 
bp, and other input quantities are taken from the 3D numerical 
simulation. When the CMP is used to obtain the other third- 
and fourth-order moments, additional quantities have to be de- 
termined, namely b and m in Eq. d46l l for Sg (see Eq. d34ll). 

5. Conclusions 

With the help of 3D numerical simulations of the upper part 
of the solar convective region, we have shown that the QNA 
and the TFM fail to describe the fourth-order velocity and 
temperature correlation moments, if merely used on their 
own. These results confirm KR2006 and geophysical studies 
( Gryanik & Hartm ann 2002) and led us to generalize the TFM 
in order to take the effects of the turbulent properties of the 
up- and downflows explicitly into account (GTFM). We point 
out that the GTFM can be used in other contexts than the solar 
one as long as the convective system can be described with two 
turbulent flows. 

One might wonder whether it is likely that the CMP and the 
model for p mode excitation developed in Paper II are generally 
applicable to solar-like stars. To answer this question requires 
further work, but results on important ingredients of these mod- 
els are encouraging. The case of convection in the planetary 
boundary layer of the atmosphere of the earth was already dis- 
cussed in GH2002. Their interpolation model for FOMs has 
mean while been inv estigated for the case of convection in the 
ocean dLoschll2004 and solar granulation dKupka & Robinson! 
20071 who also stu dy the case of a K dwarf; pr eliminary results 
were published in Kupka & Hillebrandt 2005). We corroborate 
the latter here with simulations for solar granulation based on 
more realistic boundary conditions. The overall conclusion that 
can be drawn from these studies is that, at least away from the 
boundary layers of convection zones, the FOMs in purely con- 
vective flows can be estimated according to the interpolation 
model by GH2002 with an accuracy typically in the range of 
20% to 30%, whereas the QNA is off by a factor of two to three. 
For the superadiabatic layer, the discrepancies of the QNA re- 
main the same in any case of the same size. 

We focused here on the solar case, more precisely a re- 
gion that is nearly adiabatic, just below the superadiabatic zone 
where the acoustic modes are excited. As indicated by the 
3D simulations, the coherent downdrafts, called plumes, are 
more turbulent than the upflow. In addition, we use the plume 
model developed by RZ95 to estimate the upward and down- 
ward mean velocities. With these additional approximations, 
the GTFM yields a closure model, the CMP, which can be ap- 
plied in the quasi-adiabatic zone (located just below the su- 
peradiabatic one). Comparisons of calculations based on the 
CMP with direct calculations from the 3D numerical simula- 
tions show a good agreement. Hence, the CMP provides an 
analytical closure for third- and fourth-order moments. These 
moments are expressed in a simple way and require only the 
knowledge of the second-order moments and the parameters of 
the plume model. We stress that the CMP involves four param- 
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eters: the number of plumes in the considered shell (i.e., near 
the photosphere), the exponent of the power law for the mean 
vertical velocity of plumes, the law to describe the temperature 
difference between the two flows, and the mean fractional area 
of the updrafts and hot drafts. 

A study of the dependence of the results on these parame- 
ters is in progress. For instance, an increase of a will imply an 
increase of S w in Eq. ( l3"3T l. and hence of the fourth-order mo- 
ment < w' 4 >. Nevertheless, it is extremely difficult to deduce 
the behaviour of the system, since from Eq. ( |4TT > a variation of 
a changes the velocities of the flows. Instead, one could use a 
set of numerical simulations to study the effect of a change of 
the parameter a. In a companion paper, we use the CMP in a 
semi-analytical approach to calculate the power supplied to the 
solar p modes. It is found that the power is quite significantly 
affected by the adopted closure model. 

Our final aim is to apply the CMP to the study of stochastic 
excitation of solar-like p modes in stars other than the Sun. It 
will be necessary to assess the validity of the CMP approxima- 
tions to extend their application to stellar conditions different 
from the solar case. This will also require investigating the de- 
pendence of the parameters entering the CMP, for instance, on 
the effective temperature of the star (work which is in progress). 
As pointed out in Sect. @] the CMP is valid only in the quasi- 
adiabatic zone due to the power laws used to model the plume 
dynamics. This will be discussed further in the companion pa- 
per in which the present model will be used in the superadia- 
batic zone in order to propose a new closure for the calculation 
of stellar p modes. 

Finally, we note that in the present work we do not take 
the effect of differential rotation and meridional circulation 
into accoun t . However, recent heliose i smic investigations 
( Sch ou et al.1 l2002t IZhao & Kosovichevl |2004) have shown 



Appendix A: Cross-correlation moments 



that variability of those large-scale flows gradually affects 
wavelength and frequencies, leading to a redistribution of 



the observed power spectrum (Shergela shvili & Poedt s 2005 



Hin dman et al .112005b . Hence, it could have an indirect effect 
on the amplitudes of p modes. Furthermore, large-scale 
laminar non-uniform flows can have a significant effect on the 
formation of the coherent structures and i ntrinsic turbul ence 
(Mie schet ai]|2000t iBrun & Toomrell2002t |Rempelll2005l) . To 



what extent they can affect solar p mode amplitudes, through 
the closure model and the Reynolds stresses, remains to be 
investigated. 
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As explained in Sect. 13.11 we provide the cross-correlation mo- 
ments: 

< w'ff > = a uh < w'e' >„,,, + due < w'e' > u , c 

+ a dc < w'S' >d, c + adh < w'6' > d j, + T]5w56{A.Y) 



< w' 2 6»' > = a uh < w 2 6' > uJl + Que < w a ~el >„, e 
+ a dh < w a ~& > d ,h + adc < w' 2 6' > dfi 
+ /3 l 66 + /3 2 5w + fa Sw 2 66 



< w '0' 2 > = a uh < w'e' 2 > uM + a uc < w'6' 2 >„, c 
+ a dh < w'6' 2 > dJl + a dc < w'6' 2 > dt( 
+ j\6w + jo b~e + y3 dwde 2 



(A.2) 



(A.3) 



< w e > = a u h < w 



r./2Z'2 



>u,h + a uc < W 



+ a dh < w u > d ,h + a dc < w > dyC 
+ tp\de + fa 6w + 03 66 2 + (p4 5w 2 
+ fcdwde + cf> 6 6w 2 6e 2 



(A.4) 



< VV' 2 > u 



< w 2 > d 



Pi = 2 

ft = 

J\ = 
+ 

72 = 2 

73 = 



where: 

T] - a uh (l - a)(l -b)- a uc b(l - a) + a dc ab - a dh a{\ - b) 
Pi = a uh (l -b)-a uc b 
adh(J -b)-a dc b 

a uh < w'6' >„,/, +a uc < w'e' >„ >c -a(< w'e' > -rj) 
a u h{\ — fl) 2 (l - b) — a uc b(l — a) 2 + a d i,a 2 (l — b) - a dc a 2 b 
a uh (l -a)-a dh a 
a uc (l - a)- a dc a 

a u h < w'e' >„,/, + a dh < w'e' > dJl - b(< w'e' > -rj) 
a u h{\ — a){\ — b) 2 — a uc b 2 (\ — a) — a d i,a(l — b) 2 — a dc ab 2 
h = 2 a uh < w' 2 e' >„,/, + a dh < w 2 e' > dd , 



<& 2 >h 



< e' 2 > r 



- b(<w' 2 e' >-fr - fo- fr) 



a uh < w'e' 2 >„,/, 



+ a u 



< w' 2 e' > u 



a(< w'e' 1 



7i - 7z - 73) 



93 = 



fl«A(l -b) 2 -a uc b 2 



a dh (l-b) - a dc b 



a uh (l - a) 2 -a dh a 2 



a uc (l - a) - a dc a^ 



w 2 >„ 



w 2 > d 



7*2 
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05 = 4 



a uh (l -a)(l-b)< w'& >„ ;/l - aucb(l - a) < w'& > u 



a dh a(l - b) < w'ff > d ,h + a dc ab < w'6' > d , c 



06 = 



a uh (l - a) l {\ - b) z + a uc (l - a) 2 b 2 + a dh a\\ - b) 



+ a dc a 2 b 2 
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